Analisi dei dati effettuata con la suite CONTINUUM ANALITICS https://www.continuum.io/
basata sui seguenti moduli python:
tutti i software utilizzati sono open source
In [1]:
%matplotlib inline
#importo le librerie
import pandas as pd
import os
from __future__ import print_function,division
import numpy as np
import seaborn as sns
os.environ["NLS_LANG"] = "ITALIAN_ITALY.UTF8"
In [2]:
#importo il file con i dati
path=r"D:\d\05 Lavscien\autoimmunita\corr_thibya\compar_thibya_immulite.csv"
database=pd.read_csv(path,sep=';',usecols=[1, 2, 3,4,5])#colonne da utilizzare
database['valore_cap']=database['valore_cap'].apply(lambda x: round(x,2))
database.drop(['codificato','accettazione'],axis=1,inplace=True)
database.tail(6)
Out[2]:
In [3]:
database.describe()
Out[3]:
Varibili d'ambiete in comune
In [4]:
#variabili d'ambiente comuni
cutoff_cap=2.9 #tre 2.9 r 3.3 dubbi
#cutoff_cap=3.3
#cutoff_rout=1 #brahms 1-1.5 dubbi
cutoff_rout=0.55 #Siemens
METODO_ROUTINE="Siemens Immulite 2000 Chemil."
CAP="Thermo Fisher ELIA anti-TSH-R Cap250 "
In [5]:
database['cap_PN']=(database['valore_cap']>=cutoff_cap)
database['rut_PN']=(database['valore_rut']>=cutoff_rout)
database.head(5)
Out[5]:
In [6]:
database['cap_PN'].replace([True,False],['Pos','Neg'],inplace=True)
database['rut_PN'].replace([True,False],['Pos','Neg'],inplace=True)
database.describe()
Out[6]:
In [7]:
#sci.py moduli
from scipy.stats import chi2_contingency, fisher_exact
pd.crosstab(database.cap_PN,database.rut_PN)
Out[7]:
In [8]:
ax=pd.crosstab(database.cap_PN,database.rut_PN).plot(kind='bar',stacked=True, )
ax.legend(['Neg','Pos'])
ax.set_xlabel(CAP)
Out[8]:
In [9]:
# test chi square
chi2, pvalue, dof, ex = chi2_contingency(pd.crosstab(database.cap_PN,database.rut_PN))
print ('valore di p:{}'.format(pvalue))
In [10]:
# test esatto di Fisher
oddsratio, pvalue =fisher_exact(pd.crosstab(database.cap_PN,database.rut_PN))
print ('valore di p:{}'.format(pvalue))
In [11]:
from statsmodels.sandbox.stats.runs import mcnemar
stat,p=mcnemar(pd.crosstab(database.cap_PN,database.rut_PN))
print("valore di p:{}".format(p))
-
In [12]:
# grafico di dispersione
import matplotlib.pyplot as plt
fig = plt.figure()
fig.suptitle('Scatterplot', fontsize=14, fontweight='bold')
ax = fig.add_subplot(111)
ax.set_xlabel(METODO_ROUTINE)
ax.set_ylabel(CAP)
ax.plot(database.valore_rut,database.valore_cap,'o')
plt.show()
eseguiamo ora lo studio di regressione con tre modelli diversi
Moduli statmodels e scipy
In [13]:
# con statmodel : regressione minimi quadrati
##res_ols = sm.OLS(y, statsmodels.tools.add_constant(X)).fit() per vecchia versione
import statsmodels.api as sm
#sm.OLS(Y,X)
X = sm.add_constant(database.valore_rut )
modello_minquad=sm.OLS(database.valore_cap,X)
regressione_minquad=modello_minquad.fit()
regressione_minquad.summary()
Out[13]:
In [14]:
# con statmodel : regressione robusta (Robust Linear Model)
X = sm.add_constant(database.valore_rut)
modello=sm.RLM(database.valore_cap,X)
regressione_robusta=modello.fit()
regressione_robusta.summary()
Out[14]:
In [15]:
#importo la librearia seborn per una migliore visualizzazione grafica
sns.set(color_codes=True)
ax = sns.regplot(x=database.valore_rut,y=database.valore_cap, color="g",robust=True)
ax = sns.regplot(x=database.valore_rut,y=database.valore_cap, color="b")
ax.set_title('Regressione lineare OLS + RLM ')
ax.set_xlabel(METODO_ROUTINE)
ax.set_ylabel(CAP)
ax.set(ylim=(0, None))
ax.set(xlim=(0, None))
Out[15]:
In [16]:
sns.set(color_codes=True)
ax2 = sns.regplot(x=database.valore_rut,y=database.valore_cap, color="g",robust=True)
ax2 = sns.regplot(x=database.valore_rut,y=database.valore_cap, color="b")
ax2.set_title('Regressione lineare OLS + RLM ')
ax2.set_xlabel(METODO_ROUTINE)
ax2.set_ylabel(CAP)
ax2.set(ylim=(0, 20))
ax2.set(xlim=(0, 8))
Out[16]:
In [17]:
ax=sns.jointplot(x=database.valore_rut,y=database.valore_cap, kind="reg");
ax.set_axis_labels(METODO_ROUTINE,CAP)
Out[17]:
Ortogonal Distance Regression (Deming Regression)
In [18]:
# regressione ODR (ortogonal distance regression Deming)
import scipy.odr as odr
#modello di fitting
def funzione(B,x):
return B[0]*x+B[1]
linear= odr.Model(funzione)
variabili=odr.Data(database.valore_rut,database.valore_cap)
regressione_ortogonale=odr.ODR(variabili,linear,beta0=[1., 2.])
output=regressione_ortogonale.run()
output.pprint()
output
Out[18]:
In [19]:
database_b=database
database_b['bias']=database['valore_rut']-database['valore_cap']
database_b.head()
Out[19]:
In [20]:
sns.distplot(database_b.bias)
Out[20]:
In [21]:
database.describe()
Out[21]:
Normalize data
In [22]:
'''from sklearn import preprocessing
min_max_scaler = preprocessing.MinMaxScaler()
x_scaled = min_max_scaler.fit_transform(x)
normalized = pd.DataFrame(database.valore_cap)
normalized'''
Out[22]:
In [ ]: